Calculation of reduced density matrices from correlation functions 
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It is shown that for solvable fermionic and bosonic lattice systems, the reduced density 
fS) , matrices can be determined from the properties of the correlation functions. This provides 

' the simplest way to these quantities which are used in the density-matrix renormalization 

' group method. 
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. Reduced density matrices for solvable fermionic and bosonic lattice models have been studied in recent 

years because such operators play a central role in the density-matrix renormalization group (DMRG) 
method In contrast to the quantities used in other cases, they refer to a subset of sites, not to 

a subset of particles. It has been found that they have exponential form exp(—Ti.), where Ti. is again a 
solvable fermionic or bosonic operator, confined to the chosen subsystem This was derived by 

starting from the total density matrix (usually for the ground state) and integrating out the degrees of 
freedom outside the subsystem. In the case of fermions, this can be done using Grassmann variables. 
The procedure is straightforward, but also somewhat tedious. However, it was noted recently that for 
fS| ' a hopping model the final result involves only the one-fermion correlation functions of the system 

, In the following it is shown that one can go one step further and base the considerations completely on 

correlation functions. The density matrices then follow in a very simple and transparent way. 
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' Cn.rn — < cJ^Cm > (2) 

the one-particle function in this state. The Cn,Ta form a Hermitian matrix C. Because | 5* > is a 
determinant, all the higher correlation functions can be expressed by C, e. g. 

< 4icLcfcCi > = < Clci >< c\^Ck > - < 4 Cfc >< cX,ci > (3) 

Now consider a subsystem of M sites for which the notation i,j will be used. By definition, the reduced 
density matrix p reproduces all expectation values in the subsystem. Therefore the one-particle function 

is 



Consider first a system of free fermions hopping between lattice sites. The corresponding Hamiltonian 
has the general form 



where the "hat" denotes quantities of the total system. This Hamiltonian has Slater determinants as 
eigenstates. Let | \E' > be such a state and 



= tr{pc\c,) (4) 

and the higher functions must factorize as in ). According to Wick's theorem, this property holds if p 
is the exponential of a free-fcrmion operator [p| . Thus one can write 

p = /Ccxp(-7i:) (5) 
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where /C is a normalization constant and 



Let (/)fc(i) be the eigenfunctions of H with eigenvalues e^. Then the transformation to new fermion 
operators ak 



^<Pk{i)ak (7) 



diagonalizes TL and p becomes 



M 

p = /Cexp (- ^efcaj^flfc) (8) 
fe=i 

Using this in (^) together with tr(p) = 1 gives 

C^,J ^Y.'t'li^'l'kij)^^ (9) 
k 

On the other hand, H has the representation 

H^^, ^Y.^ki^)m) (10) 

A: 

Therefore the eigenvalues of the two matrices are related by 

a = {e^" + 1)-' (11) 

and in matrix form, with the prime denoting the transpose 

H' ^ ln[{\ - C)/G] (12) 

This is the formula found in Q. Due to its form, p is completely determined by the M x M matrix 
C. One should note that any one-particle correlation function can be expressed in such a way through 
a proper free-fermion operator. The only condition is that the eigenvalues Cfc of C lie between and 1 
and this is always the case, since they can be written in the form < a^Ofc > with new fermion operators 
However, for a state which is not a Slater determinant, the free-fermion density matrix found above 
would in general give wrong results for other expectation values. 

These considerations can be extended to systems with pair creation and annihilation which can be 
diagonalized by a Bogoliubov transformation. The eigenstates are then Slater determinants in the new 
Fermi operators. In such a state, "anomalous" correlation functions 

Fn.m —< '^n'^ln (13) 

exist which also occur in the factorization equations. Thus (^) is changed into 

< C^cl^CkCl > = < C^Cl >< 4jCfc > - < 4cfc >< cJ„Q > -t- < >< CfcQ > . (14) 

To reproduce this, p has to be an exponential with an operator 7i which also contains pair creation and 
annihilation processes 
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W = Et^I^^^-^^- + + h.c.)] (15) 



Since now two matrices appear in Ti, one needs additional input, which is provided by the correlation 
functions Fi j. By following the usual diagonalization procedure for H, [0 and calculating C and F one 
can then show that 

[{C - 1/2 -F){C~ 1/2 + F%,, = \Y. M^WkU)th^ iek/2) (16) 

k 

where the (j)k{i) are the orthonormal eigenfunctions of {A — B){A + B) and the £k are again the single- 
particle eigenvalues of TL. Thus one can find the Efc from the eigenvalues of the matrix on the left hand 
side of (|l6|). This matrix can be written as KK^ /A where K/2 = (C— 1/2 — F), since F is anti-Hermitian. 
For _F = 0, the result is equivalent to (^. If one turns a hopping model into one with pair terms via a 
particle- hole transformation, (^) goes over into (|l^). To make contact with the treatment in |^, one first 
relates the matrix G, used there to write an even-parity eigenstate in the form 

|vI/>=Cexp{i^G„„44}|0>, (17) 



where | > is the vacuum of the c„, to the quantities appearing here. For the ground state, G connects 
the two sets of functions t^p and "ipp arising in the diagonalization of the Hamiltonian. The same holds 
for k/2 = {C -112- F) and one finds that G = [K - l)l{k + 1). Using this, one can show that the 
matrix in (A9) of Q equals 2(1 + K K^) / { \ — KK'^) jll]. Therefore the eigenvalue equation used in |^ to 
determine 2chek is an alternative version of the relation ( p^ and both approaches are fully consistent. 

In a similar way, one can treat systems of coupled harmonic oscillators. In this case, it is convenient 
to consider the correlation functions of positions and momenta 

In the ground state, which is a Gaussian in the coordinates, one then has factorization formulae like 

< XnXmXkXl > = < XnXm >< XkXl > + < XnXk >< XmXl > -f < XnXl >< XmXk > (19) 

which are non-trivial even if all indices are equal. They hold also if the expectation values are calculated 
with an exponential operator quadratic in the x and p. Therefore p has the form (^) with 

7^ = i E [T^jpipj + Vi,jXiXj] (20) 

The diagonal form is again (||) but with bosonic operators and the ek follow from the eigenvalues i/^ of 
the matrix XP via 

cth{ek/2) = i^k/2 (21) 

If the subsystem is a single oscillator i, there is only one ly given by < a;f >< > . For a homogeneous 
system, this can also be expressed through the frequency moments of the normal modes as < l/o; >< >. 
The general equations are again equivalent to those obtained previously in ||^. 
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This shows that the way to reduced density matrices associated with eigenstates of solvable fcrmionic 
or bosonic systems can be shortened considerably. The results are also valid for systems at finite tem- 
perature as considered in 0|. In connection with the DMRG, the main aim has been to determine the 
spectra of the p and their general features. Here the present approach helps if the necessary correlation 
functions have simple analytic expressions. This is the case for nearest-neighbour hopping on a chain, 
or on a square lattice for half filling. Then only the diagonalization of the matrix C remains. There are 
some limitations, because if large Sk occur, the corresponding eigenvalues of C are exponentially close 
to or 1. Also the relation of p to the corner transfer matrices of two-dimensional models is not 

visible unless one determines TL in its non-diagonal form. Still, one comes very close to an analytical 
solution, and working with the correlations can give additional insight into the nature of the problem. 

Note added : Formulae similar to those given here can also be found in a recent paper on entangled 
quantum states fl^ . 
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